Chemical Oscillations out of Chemical 

Noise 



Carlos Escudero 

Departamento de Economia Cuantitativa 
<fe Instituto de Ciencias Matematicas (CSIC-UAM-UC3M-UCM) 
Universidad Autonoma de Madrid, 28049 Madrid, Spain 
e-mail: cel@icmat.es 

Andres M. Rivera 

Departamento de Ciencias Naturales y Matematicas, Facultad de Ingenieria 
PontiGcia Universidad Javeriana Cali, 26239 Cali, Colombia 
e-mail: amrivera@puj.edu. co 

Pedro J. Torres 

Departamento de Matematica Aplicada, Facultad de Ciencias 
Universidad de Granada, 18071 Granada, Spain 
e-mail: ptorres@ugr.es 

Abstract 

The dynamics of one species chemical kinetics is studied. Chemical 
reactions are modelled by means of continuous time Markov processes 
whose probability distribution obeys a suitable master equation. A 
large deviation theory is formally introduced, which allows develop- 
ing a Hamiltonian dynamical system able to describe the system dy- 
namics. Using this technique we are able to show that the intrinsic 
fluctuations, originated in the discrete character of the reagents, may 
sustain oscillations and chaotic trajectories which are impossible when 
these fluctuations are disregarded. An important point is that oscil- 
lations and chaos appear in systems whose mean-field dynamics has 
too low a dimensionality for showing such a behavior. In this sense 
these phenomena are purely induced by noise, which does not limit 



itself to shifting a bifurcation threshold. On the other hand, they 
are large deviations of a short transient nature which typically only 
appear after long waiting times. We also discuss the implications of 
our results in understanding extinction events in population dynamics 
models expressed by means of stoichiometric relations. 



1 Introduction 

The kinetics of reaction systems have been extensively studied over the years. 
These systems have been postulated as paradigmatic models for the de- 
scription of a large number of natural phenomena, including topics from 
organic and inorganic chemistry, epidemiology, population biology and gene- 
tics, nuclear physics, non-equilibrium statistical mechanics and many others 
sciences [HE]- The mathematical description of such systems usually starts 
with the assumption of a set of stoichiometric relations of the form 

A^B, (1) 

signifying that the reagent A transforms to B with the time dependent proba- 
bility 

PA^B{t) = ae-"*, (2) 

where a > is the specific reaction rate. Note that the probabilistic nature of 
the reactions introduces fluctuations into the dynamics: this is the "chemical 
noise" we will be interested in. In more general terms, the state of a system 
containing m reagents and n reactions is described by a continuous time 
Markov process. All the available information is carried by the distribution 
P{nA,nB, ■ ■ ■ ]t) specifying the probability of the existence of exactly ua 
reagents of type A, ub of type B, at time t. The dynamics of this 
distribution is governed by a master equation of the form [H [2] 

dP 

where Pk denotes the probability of finding the system in the state k charac- 
terized by a certain number of reagents of each type, and Wj^k is the tran- 
sition matrix from state j to state k. While solving the master equation 
to find Pk would mean that we possess all the available information on the 
system, the chances of obtaining an exact closed form for Pk are scarce in 
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realistic situations. Additionally, the amount of information it brings is usu- 
ally excessive, and a great part of it does not add any valuable information 
about the dynamics. Consequently, the most common approach to the sub- 
ject concentrates on the dynamics of some statistical quantity of interest, 
as for instance a density, which is able to describe the system macroscopic 
state. The selection of an adequate variable has to be supplemented with 
selecting a suitable approximation in order to get an operative theory that 
allows studying the otherwise commonly untractable master equation. A 
particularly advantageous choice is the analog of the quantum mechanical 
Wentzel-Kramers-Brillouin (WKB) approximation adapted to this sort of 
systems, which is now well established in both physical and mathematical 
literatures, see for instance [3llilEll6l[7llHll9l[Tni[IIlll2l[13l[Tll[T5]. It allows 
the description of both the short time dynamics, which is to a large extent 
independent of the fluctuations and therefore captured by mean-field type 
approximations, and the long time behavior which is affected, dramatically 
on occasion, by large deviations. The mathematical and physical natures of 
these large deviations will appear clearer in the following sections. 

In this work we are concerned with simple reaction sets which on the 
other hand have an intuitive physical meaning. This way we will focus on 
simplified model systems which, although not of direct practical applicability, 
facilitate analytic progress and physical intuition. We will show how chemical 
fluctuations strongly affect the dynamics for long times, when large deviations 
have had time to develop. In these cases, chemical fluctuations are able 
to promote periodic orbits and chaotic behavior in reaction systems whose 
dimensionality is too low to present such a behavior if strictly deterministic 
dynamics are considered. These effects are, however, both rare and short- 
lived. They are, at the same time, substantially different from other sorts 
of noise-induced oscillations which appear in different important phenomena 
and whose structure relies on an existing deterministic mechanism (like the 
proximity to a deterministic bifurcation) which is anticipated or activated 
by noise [121 [IZl [THl [IHl 120] • In this sense, we may say our focus is on 
oscillations which are purely promoted by chemical noise. Our approach will 
be probabilistic at the beginning, when we will present formal calculations 
in which the equations governing large fluctuations will be derived. These 
equations have the form of Hamiltonian dynamical systems, which will be 
in realistic situations genuinely different from the ones usually considered in 
classical mechanics. For them we will be able to show, by means of explicit 
calculations, rigorous proofs and numerical simulations, how chemical noise 
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is capable of sustaining chemical oscillations, both of periodic and chaotic 
nature. 



2 Large deviations 
2.1 Brownian motion 

We devote this section to clarify the type of problems we are going to solve. 
We start with the perhaps simplest stochastic process one could consider: 
one-dimensional Brownian motion. For our current purposes it will be the 
solution of the equation 

f-om, (4) 

where D > is a diffusion constant and C,{t) is the standard Gaussian white 
noise defined by its two first moments 

where 6{-) is the Dirac delta distribution. Of course, a rigorous interpretation 
of this equation is possible in terms of Ito calculus [21], but such a precise 
definition will not be needed herein. Equation (|4]) is provided with the initial 
condition B{0) = 0. A classical problem within this subject is calculating the 
first time the random walker B{t) reaches some fixed level a 7^ in absence 
of other constraints. The well-known solution states that the random walker 
reaches level a in finite time with probability one, but the mean time at which 
this event occurs diverges. 

Langevin equations like ^ and more complicated variants of it are well 
understood from the large deviations point of view [22] • It associates to this 
stochastic differential equation the rate or action functional 

SHt)] = ,^ [ \x{t')\'dt', (6) 



which in the small noise limit D — t- gives rise to the following Euler- 
Lagrange equation 

X = 0, (7) 

for the position x of the random walker. If we complement this equation 
with the boundary conditions x{0) = and x(T) = a we find the solution 

= (8) 
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signaling the most probable trajectory that links the origin with the level a 
after a time T for the Brownian dynamics Q. Small deviations are obtained 
by setting D = in Q. In this case the random walker stays at the origin 
for all times. So the full picture would be, for small noise, the random walker 
will be at a neighborhood of the origin with a large probability but with a 
small probability large deviations might appear and drive the system further 
away. The probability V with which these large deviations, which promote 
trajectories ([s]), manifest themselves into the system dynamics is proportional 
to the exponential of the negative of the action 

P~e- = exp(-^). (9) 

From this formula it is clear that those trajectories that connect the origin 
with the level a in a shorter time are rarer than those which take a longer 
time. The prefactor in this case is easily found by normalization. Note that 
the large deviation theory is valid for ^ DT, otherwise the system diffuses 
away from the original position and the approximation breaks down. 



2.2 Plankton extinction 

We will describe now the large deviations technique in the context of reac- 
tion kinetics. To this end we consider a simple model that has nevertheless 
a genuine practical interest. This model was introduced as a description 
of plankton population dynamics |23l EH [25] and nonequilibrium statistical 
mechanics p6j. It consists of the following reactions 

A^2A, A ^ 0, (10) 



happening at the same rate 7. We will employ large deviation theory to 
describe the extinction of the "plankton population" A. The state of the 
system may be described by a continuous time Markov process obeying the 
master equation 

dP 

= 7[(n - l)P„_i - nPn] + i[{n + l)P„+i - nPnl (11) 

where the term inside the first bracket corresponds to the branching reaction 
and the one inside the second bracket to the disintegration reaction. By 
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introducing the generating function 



G(p,t) = 5^p"P„(t), (12) 



n=0 

we transform the master equation into the following "imaginary time Schrodinger 
equation" 

dtG = ^ip-lYdpG. (13) 

Note that we are employing the "momentum" rather than the "coordinate" 
representation in this last equation. One can obtain the probability distri- 
bution from the generating function in the following way 

Pn{t) = ^jGip,t)p--^, (14) 

where the contour integral runs over a closed path on the complex p— plane, 
surrounding the origin and going through the region of analyticity of G{p, t). 
The corresponding "classical" Hamiltonian of our theory reads 

n = i{p-lfq- (15) 

it is precisely this Hamiltonian, as in the previous case, which describes the 
large deviations of the system. The corresponding equations of motion are 

q = _ = 27(p-l)g, (16) 
rH-l 

p = -^ = -l{p-^f- (17) 

Note that the line p = 1 is degenerate and all points on it are fixed points. 
These solutions refer to small deviations: the system stays with a large prob- 
ability in a neighborhood of the initial condition for short times as in the 
Brownian motion case. The solution for the coordinate q is 

q{t) = [vW)±ViHty, (18) 

where the minus sign is selected for extinction trajectories and if is a con- 
stant indicating the initial "energy" . The duplicity of signs in this equation 



comes from the extraction of the square root of equation (15). The number 
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n of reagents can be calculated by means of formally applying a steepest 
descent approximation to formula (14) for G{p,t) = exp[—S{p,t)], where S 
is the "classical" action [9]. Then one finds n ^ —pdpS, and employing the 
"classical" relation q = —dpS one concludes n{t) ^ p{t)q{t). In our particular 
example one finds 



n 



pq = q 




(19) 



which becomes zero due to a fluctuation when p — )■ 0, which leads to g = H/ 7. 
If the system follows an optimal trajectory it will become extinct after a time 



lq{0) 



jH 7 



(20) 



The probability with which this realization appears for short times is the 
exponentiated negative of the action V ~ , up to some prefactor. We will 
limit ourselves to the exponential order, as the calculation of the prefactor is 
a rather technical issue [15] and will not add substantial information to the 
present discussion. The action reads 



S 



(pq - H)dt + [p(0)g(0) - p{tMte)] + S^ = So, 



(21) 



in the case of an extinction trajectory, where the initial action 5*0 = — \n[G{p, 0)]. 
The last equality has been derived employing the following derivations 



{pq - H)dt 



d 
di 



(pq) 



[pq--f{p- lfq]dt 



7g(/-l)rft, (22) 
(23) 



where we have substituted p and q for their respective values from (16) and 



(17). We consider two initial conditions as in [9], the Poisson distributed 
initial condition with average no, this is (7(^,0) = exp[no(p — 1)], and the 



Kronecker delta centered at ng, which is G{p, 0) 
extinction probability reads 



p^o. In the first case the 



V ~ exp 



H"^ Hno H 

+ + 



47 



7 



27 



(24) 



7 



and in the second 



V 



H 



+ 



H 



no 



and both yield the same result in the thermodynamic limit no — ?• cxo 

V ~ exp 




(25) 



(26) 



The optimal trajectory corresponding to this characteristic time to extinction 



is found by combining Eq. (18) (with the minus sign) and the expression for 



n given by the second equality of Eq. ( 19 ) 

n{t) = no + 'jHf - t^/lP^^An^ 



(27) 



where the thermodynamic limit no — oo has been considered in the last 
step. In the derivation of the first equality we have employed the following 
two relations 



nit) 



1 - 



W) 



t 



g(0) 



no 



+ 4 no. 



(28) 



(29) 



V ^ 

Note that we have found a one parameter family of solutions, parameterized 
with the energy H. Time is not the mean extinction time, because at ev- 
ery time there are equally probable trajectories which do not become extinct. 



the ones corresponding to the plus sign in Eq. (18). This makes the mean 



extinction time infinite, although the system becomes extinct with proba- 
bility one [21 [27] . The interpretation of this time is that of a characteristic 
time of extinction, this is, if the system becomes extinct after a time then 



the most probable path to extinction would have been (27). Note that more 



"energetic" trajectories lead to extinction faster but they are less probable. 
Using the relation between te and H we may find an expression akin to ([9]): 



V ~ exp ( — 



no_ 



(30) 



Note that in this case the scaling is different. Note also that, as in the 
previous case, the large deviation theory that has led us to the optimal 



trajectories (27) is valid for short times t ^ no/7. 
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3 Chemical oscillations and chaos 
3.1 Noise induced oscillations 



T ^ ^ ^^^^^^^ ^ V i 

f / V V J 

^ , , ; 

1 1 1 ^ - - ^ ^ M t 

, ' ^ ' ^ 

; , X ^ -^-^-.^ ^ ^ ^ + 

r K ^ ^ ^ i 



Figure 1: Vector field of tlie Hamiltonian system (36) (detail of tlie qua- 
drangular area enclosed by the four zero energy lines). The values of the 
parameters are = 1 and a = 2. 



We now move to studying the more complex nonlinear reversible reaction 

A 2A, 2A ^ A. (31) 

It is clear that it can be considered as a stochastic discrete model for logistic 
growth. The master equation describing this reactions set is 

^ = /i[(n - l)P„_i - nPn] + ^[in + l)nPn+i - n{n - 1)P„]. (32) 

We may use the generating function technique to convert this equation into 
a partial differential equation which can be exactly mapped, using quantum 
mechanical tools, into the path integral [9] 

U = [ PpPge-^IP'"!, (33) 
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where U is the problem Green function. The action reads 



pq - - p)g - 2 (P - P^)^^ + [p(0)g(0) - p{t)q{t)] + Sq, 

(34) 

and upon rendering the variables nondimensional q — > {2fi/a)q (so this new 
q = 0(1), see below) and t ^ t/fi one finds 

2u f Z"*/^ 

^b, ^] = ^ W N - (P' - ^^)^ - - P')?'] + [^^(O)g(O) - p{t)q{t)] + So 

(35) 

where 5*0 = aSo/{2fi) = 0(1), so the steepest descent method makes sense 
for a. In this approximation and recovering the dimensional coordinates 
the large deviations problem reduces to studying the Hamiltonian ^llj 

'H = fiip"^ - p)q + ^{p ~ p'^)q^ , (36) 



and the corresponding dynamical system 

on 



dq 



— = li{2p-l)q+ -{l-2p)q\ 



(37) 



This system has five fixed points, four of which lie in the energy H = level, 
these are 

(0,0), (0,2/i/a), (1,0), and (1,2/i/a), (38) 

all of them are saddles. The fifth fixed point is {1/2, fi/a), its energy is 
H = — /i^/(8(T), and it is a local minimum of energy, what implies it is a 
center. The H = level is particularly simple, as it is composed of the four 
invariant lines 

{p = 0}, {p=l}, {q = 0}, and {q = 2fi/a}. (39) 

The dynamics is exactly integrable in all these lines. They cross at the four 
zero energy fixed points, and they enclose a quadrangular area in whose 
center lies the fifth fixed point. In this quadrangular area all the trajectories 
are periodic orbits surrounding the center, see Fig. [TJ Note that, both outside 
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and inside these four zero energy lines the energy is strictly negative. The 
mean-field behavior corresponds to the {p = 1} line, on which the dynamics 
reduces to the well known logistic equation 



(i = m- (40) 



As expected, reaction (31) corresponds to pure logistic growth when the 
fluctuations are neglected, this is, for a large number of reagents and short 
times. So at the mean-field level the only possibility is a monotonic approach 
to the stable fixed point q = 2fi/a, provided the initial condition fulfills 
g > 0. So the dynamical scenario presents a very reduced phenomenology in 
this case. However, if we consider the intrinsic fluctuations and thus the full 
phase space things are different. In this case, for instance, we may observe 
the system for a short time in the neighborhood of the flfth flxed point, which 
represents a reagent density n ~ /i/(2(T) (this is obtained as the product pq 
evaluated at the flxed point). If we initialize the system with this reagent 
number we have a probability 



T> -s[p,q] / \ ^ r /^^^ /iln(2) 
V ^ e = exp , V ^ exp 



(41) 



of observing the system in a neighborhood of this point during a time r 
respectively for the Poissonian distributed and deterministic initial condition. 
But furthermore we can observe periodic behavior. All the periodic orbits 
are optimal trajectories that can be observed experimentally if we wait long 
enough. These orbits are characterized by an energy > Hp > — yU^/(8a), 
and therefore the probability of observing a number m of cycles is 

again at exponential order, where tp is the time it takes to perform one such 
cycle and A is the phase space area enclosed by one such trajectory. We 
expect that formulas like this will be able to express the order of magnitude 
of the corresponding probability, not just an exponential dependence, as we 
have observed in other cases when the system is not in the proximity of 
an absorbing state |1^. So we see that, while the mean-fleld description 
predicts monotonic approach to a stable flxed point, the stochastic theory 
allows the appearance of transient periodic trajectories. Let us emphasize 
that such trajectories are not the typical behavior of the solution to the 



11 



master equation. They are large deviations, which manifest themselves only 
after very long times and are of a short transient nature. These characteristics 



are quantitatively described by the small probability of its occurrence (42) 



Of course, periodic orbits in the (p, g)— plane are not of physical nature. But 
it is on the other hand immediate that the number of reagents n{t) = p{t)q{t) 
is periodic if both p{t) and q{t) are periodic. We have represented the time 
evolution of n for an initial condition belonging to one of the periodic orbits 
in the {p, g)— plane in Fig. 12} 




Figure 2: Number of reagents n{t) = p{t)q{t) versus time t obtained from 
numerically integrating dynamical system (37). The system is initialized with 
the conditions p(0) = 1/4 and g(0) = 3/4 which correspond to n{0) = 3/16. 
The values of the parameters are n = 1 and a = 2. 



The fact that Hamiltonians like (36), which come from a chemical master 
equation, are not hermitian translates into the impossibility of expressing 
probabilities like (42) in terms of the physical variable n rather than the 
formal auxiliary variables p and q. Despite this undesirable fact, we can still 
characterize periodic orbits like the one represented in figure |2] by means of 
its period. Indeed, it is clear that for periodic solutions, the period of p{t) 
and q{t) will uniquely determine the period of n{t). So a way to connect the 
physically measurable quantity n with formula (42) is through the period of 
the oscillations of n. Of course, together with these large deviations, small 
fluctuations will be present all of the time. A way of distinguishing both of 
them is by means of their amplitude. The amplitude of small fluctuations is 

O (^\/fi/(J^ while the amplitude of these oscillations promoted by large devia- 
tions is O (/i/cr). So the difference should certainly be measurable in the limit 
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jj, ^ a, which is exactly the range of vahdity of our WKB approximation. 



3.2 Chemical chaos 



Not only periodicity but also chaotic trajectories are possible in this simple 
system. To obtain them we allow the branching rate to be a periodic function 



of time /i — > > 0. Note that, due to the structure of system Eqs. (37), 
we could let either /x, a or both be time dependent and still reduce the system 
to a /X time dependent one (while a remains constant) by means of a change 
of the temporal variable. In this case we deal with the system 



P 



fi{t){p-p'^) + (r{p^ -p)q, 



a 



/i(t)(2p-l)g + -fl 



2p)q' 



(43) 
(44) 



Herein we still can identify three invariant lines: {p = 0}, {p = 1} and 
{q = 0}. The mean-field dynamics appears on the {p = 1} line, and is 
expressed by the equation 



q = fi{t)q 



a 



r9 



(45) 



This differential equation is of Bernoulli type and can be solved to yield 



<l{t) 



g(0) exp 






l + |g(0)£exp 




dti 



(46) 



Assuming that /i(t) = /i + e/i(t), where h{t) is T— periodic and continuous 
and e is a small enough constant, we know that there exists one T— periodic 
solution qs{t) which attracts all initial conditions g(0) > 0. This is the 
solution whose initial condition fulfills 



g(0) 



exp 


Jq fi{ti)dti 




1 


f /o^exp 




dti 



(47) 



The points (0, 0) and (1, 0) continue to be fixed points in the non-autonomous 
system, and the point {0,2fi/a) gives rise to a periodic orbit on {p = 0}, 
which is formally identical to qs{t), but it is unstable on this line, and we 
will refer to it as qu{t). While the invariant line {q = 2fi/a} is not present 
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in the perturbed system, the periodic trajectories qs(t) and quit) are still 
connected. They correspond to fixed points and of the Poincare map 
associated to the forced continuous dynamical system. One may apply the 
theory developed in [27] to see that for generic perturbations h(t) the unstable 
manifold of Qs intersects the stable manifold of Qu and thus guarantees the 
existence of a heteroclinic connection linking both fixed points. 

The irregular behavior of the system dynamics comes from the fact that 
the periodic trajectories in the autonomous system may become quasiperiodic 
or even chaotic in the perturbed one. This can be justified by classical 
arguments from the Kolmogorov-Arnold-Moser (KAM) theory flEl [22] as 
follows. 

Let us consider the stable equilibrium P = (1/2, fi/a). The Floquet multi- 
pliers are by definition the eigenvalues of the corresponding Poincare matrix. 
By the Hamiltonian structure, the Floquet multipliers are complex conju- 
gate numbers Ai, A2 such that A1A2 = 1. For P, a direct computation on the 
linearized system gives Ai = e^^'^, A2 = Ai, with u = fi/2. The equilibrium P 
is said to be non-degenerate if uT 7^ kn, for k = 1, . . . ,4. A non-degenerate 
equilibrium is persistent under small perturbations as a fixed point of the 
Poincare map. In other words, P is continued as a T-periodic solution of the 
perturbed system (43)-(44) for small values of e. Besides, the presence of the 
heteroclinic loop corresponding to the energy level H = in the unperturbed 
system guarantees that the center around P is not isochronous, that is, the 
Poincare map is of twist type. Under such conditions, a generic perturbation 
gives rise to a classical KAM scenario, composed by a dense set of invariant 
tori (corresponding to quasiperiodic solutions) that are gradually destroyed 
as the perturbation increases, giving rise to domains of chaotic motion (Smale 
horseshoes) intermingled with stability islands. 

Figure [3] shows a chaotic orbit surrounding a set of five stability islands. 
Such chaotic orbits arises from the destruction of an invariant torus that per- 
sists until a critical value of the perturbation parameter e. Figure |4] presents 
a zoom of the latter picture, where the typical fractal structure can be ap- 
preciated. 

Let us mention that the stability islands are centered in periodic orbits 
of higher periods (or subharmonic solutions). In the case of Figure [sj a 
subharmonic solution of order 5 is located at the stability islands. Section 5 
is devoted to the study of the existence of such subharmonic solutions. 

We finish this section saying that we can compute the probability with 
which a chaotic orbit appears by means of the exponentiated negative of the 
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Figure 3: Numerical simulation showing chaos for system (43)- (44) with 
h{t) = sin {2t),fi = 1, 0" = 2, e = 0.97. It is drawn the Poincare section of a 



single orbit with initial values p{0) 
have been computed. 



0.53, g(0) = 0.91. More than 10^ points 



action. This has already been done in explicitly time dependent chemical 
systems for the simpler extinction trajectories |30]. Herein we have shown 
that for exponentially long times the sort of Hamiltonian chaos we have de- 



scribed is possible for the simple reaction (31) by virtue of chemical noise. 
In absence of noise only periodic trajectories are possible. We have rep- 
resented in figure |5] the resulting aperiodic trajectories for the number of 
reagents n{t) = p{t)q{t) for two different time slots and initial conditions. 
Their physical interpretation is analogous to that of the periodic case in the 
last section. It is interesting to note that some of the irregular motions that 
can be observed in certain stochastic reaction dynamics might have an un- 
derlying deterministic structure (we are always referring to large deviations); 
let us recall that the time evolutions represented in figure [5] are purely deter- 
ministic. 

4 Global continuation of the equilibrium point 



In subsection 3.2 we point out that the equilibrium point P = {1/2, ^/cr) of 
system (37) is not degenerate if uT ^ kn, for k = 1, ... ,4, with u = /i/2. 
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Figure 4: Zoom of Fig. |3] Besides the big stability islands, some smaller 
holes can be appreciated, conforming to a fractal structure. 



This property implies the continuation of P as a T-periodic solution of the 



perturbed system M3|-(44) for small values of e. In this section we find an 



explicit interval [0, E*] for the parameter e where this continuation exists and 
is unique. The main result is inspired by [31], where a similar technique is 
applied to the classical pendulum equation. 

For simplicity in the calculations we apply the following change of variables 

(48) 



to the perturbed system (|43|)-(|44|) 
T : 



-> 11 

and obtain the new system 

P = - 
4 



Tip, q) 



P 



1 

Z cr 



■^q + (yp^q + eh{t){^ 



P 



—p — aq^p + 2eh(t){q + —)p 
a a 



(49) 



Note that for the perturbed system (49), the invariant lines are 
{p=-l/2}, {p=l/2}, {g = -^/a}. 

In the following, we assume that uoT /2 ^ N with u = fi/2. 

Theorem 1. . There exists (5 = P{uj,a) such that for e G [0,-E*[ with 
-, with h* = maxtg[o,T] {h(t)}, system (49) has a unique non- 



E* 



a 



iujf3h* 
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10 20 30 40 50 
(a) 





(d) 



numerically integrating dynamical system (|43|)-(|44|) 

2 and e 



Figure 5: Number of reagents n{t) = p{t)q{t) versus time t obtained from 

The values of the pa- 
sin (2t). Panels 
=1/5 and 
The system 



and (b) 



rameters are = 1, a = 2 and e = 0.97; the function h(t) 

The system is initialized with the conditions p{0) 

2/25. Panels 
3/5 and g(0) 



g(0) = 2/5 which correspond to n{0) 
is initialized with the conditions p{0) 
to n(0) = 3/20. 



and (d) 



1/4 which correspond 



17 



trivial T-periodic solution $(t,e) = {(j){t,e),ip{t,e)) as a continuation of the 
equilibrium point P = (0, 0) of the autonomous case. 

Remark. From the proof, f3 = f3{u, a) is explicitly given by 



l3 = l3(uj,a)= max / \G(t,s)\ds, 
te[o,T] , 



where G is a Green's matrix function associated to the system (49) defined 
below, I . I means the usual uniform matrix norm. The explicit bound 

T r 1 2 a 



is easily derived. 



^(w, ^) < ^^^^T77^ max , 

I sm[ul / 2j\ [2 a 8uj 



Proof. We rewrite the system ( 49 ) in the form 



X = AX + B{t,X,e) 



(50) 



with 



and 



A 



B,{t,X,e) = cxp\ + eh{t){^-p'), 

B2{t,X,e) = -aq^p + 2eh{t){q+^)p. 

a 



Now if X{t,e) is a T-periodic solution of (50) the method of variation of 
constants provied us the next formula 



X{t,e) 



G{t,s)B{s,X{s,e),e)ds, 



(51) 



where G{t, s) is the Green's matrix function associated to this problem given 
by 

fj-V*-^)^; z/ 0<s<t<T 
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where the matrix J is defined by 



1 — cos uT 



sin uT 



sin uT 1 — cos uT 



a 



Exphcitly, the Green's matrix is 

/ sinu;(T/2 - (t - s)) 



o cosu;(T/2 - {t - s))\ 



sin(wT/2) 
4a; cosu;(T/2 - {t - s)) 



Au sm{uT/2) 
smuj{T/2 - (t - s)) 



\ a sm{ujT/2) 
for all < s < t < T and 

/ sincc;((s - t) - T/2) 

G{t,s) = l 



sm{uT/2) 
Aco cosuj{{s-t) -T/2) 



sin(wT/2) 
a cosu}{{s-t) -T/2)\ 

I 



Au sm{uT/2) 
sinw((s-t) - T/2) 



sin(u;T/2) 



\ a sin(wT/2) 
for < t < s < T. In Appendix 1 we present the explicit calculation of (51 ). 

Consider = M x [0,i?*] and the normed space 

E = {X e C{n,R'^) : X IS T -periodic] , 
with the norm \\-\\^- We define the operator T : E ^ E given by 

{TX){t):= [ Git,s)B{s,X{s,e),e)ds, 
Jo 

which is a completely continuous operator (with the norm \\-\\^ ) from E to 
itself. It follows from (51) that X{t,e) is a T-periodic solution of (50) if and 
only if X is a fixed point of T. 

Now we concentrate on estimating a value E* where the operator T will be a 
contraction a let invariant a closed ball B = {X G E : < p} for some 

positive number p = p^u, a, e). To this end, let X = {p,q), Y = {p, q) inside 
the ball E. Observe that 



\\{rx){t) - [TYmi 



G{t,s){B{s,X,e) - B{s,Y,e)ds 



< f max [ \G{t,s)\ds]\\B{s,X,e) - B{s,Y,e) 
\te[o,T] Jq J 
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where 

B — Bill, a) — max / \Git,s)\ds— max f max / \Gi At, s)\ds] . 

Next we consider 

B,{t, X, e) - Y, e) = e/i(t)(p^ - p') + a{p'q - fq). 
Notes that 

\eh(t)(p^ - p'^)\ < 2eh*p \\X - Y\\^ , with h* = max {h{t)} , 

and 

Ip^Q-P^qI = \{p-p){pq+pq) +P^{q-q)\ < 3p^ ||X - 
therefore 

\\B,{t, X, e) - B,(t, Y, e) \\^ < 2p{eh* + 3ap/2) \\X - Y\\^ . 
On the other hand 

\B2it,X,e)-B,{t,Y,e)\ = \2ehit)iqp-qp} + a{fp- q'p) + {p-p)\. 

a 

For the first two terms in the right hand we have that 

\q^p-q'^p\ = \{q- q){q+p)p + q^p-pq\ < 3p'^ \\X - Y\\^ , 

and 

\qp -qp\ = \p{q - q) + {p - p)q\ < 2p\\X - Y\\^ 
In consequence 

\\B,{t, X, e) - B2{t, F, e) \\^ < (2p{eh* + 3ap/2) + 2eh*{p + p/a)) \\X - Y\\^ 
This estimate imphes that 

\\B{t, X, e) - B{t, r, 6) 11^ < (2p{eh* + 3ap/2) + 2eh*{p + p/a)) \\X - Y\\^ 
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Finally 

\\TX - TY\\^ < p(2p{eh* + 3ap/2) + 2eh*{p + p/a)) \\X - Y\\^ 

Prom this last inequality, we concluded that the operator T will be a con- 
traction if ^ 2 1 

-ap^ + Aeh*p + -eh*p < -. 
2 a p 

Now we use the fact that the quadratic convex function 

f{p) = I'^P'' + ^^h*p + -eh*p - \, 
Z a p 



is negative in [0, po[ if /(O) < 0, (this is equivalent to have e < - — — — ) with 



Ap^h* 



Po 



3(7 



In this way for e e [0, E*[ with E* 



we have 



\\rx-rY\L<k\\x-Y\\^ 

with < /c < 1 for all X, F e B. Remains to knows the size of radius of the 
ball B. Let Yq = (0,0), then 



G\ i{t, s)h{s)ds 



G2,i{t, s)h{s)ds 



moreover 



lirxii^ - iirrolloo < \\rx - n\\^ <k\\x- Yo\ 



this implies 



\\rX\\^ < k \\X\\^ + IITTolL -kp + eh* (3/4. 
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eh* 13 

We take p = — ^ and e small enough sucht that /(p) < 0. Combining 
these estimatives we have that 

\\rX\\^<p, and \\TX-TY\\^<k\\X-Y\\^, 

for all X,Y E B. Then T maps the closed ball B into itself. Thus it follows 
from the Schauder fixed point theorem [33] that T has a fixed point $ in B. 
Since T : B ^ B is a. contraction then T the fixed point is unique. □ 



5 Local continuation of periodic solutions from 
the autonomous case 

From Section 3 we know that P = (1/2,/i/cr) is an equilibrium point for 



the autonomous system (37). Moreover P is a center and therefore there 
is a domain of periodic solutions around P. In Section 4, after a change 
of coordinates ( [48|, t he point P is extended for e GjOji?*] as a T-periodic 



solution e) of (49 ). Furthermore, this extension is unique and continuous. 



In this section, we look for multiplicity of periodic solutions. To this purpose. 



we assume that h{t) is an even function. Under this assumption, (49) has 
the following symmetry 

{t,P,q) {-t,-p,q). 

Our aim is to obtain nT-periodic solutions of (49) as a result of the local 
continuation of nT-periodic solutions of the autonomous Hamiltonian system 
(37). To this end we present the next result which is inspired by the results 
on [32i Section 5]. 

Some notation is needed. [■] denotes the integer part function. Fix n^, = 
r ]. For a fixed integer n>n^, define "dn = [" ^ ^ ]. 

Theorem 2. For all n > n^, and m = 1, . . . there exists en,m > such 
that for all < e < en,m system has a non-trivial nT-periodic solution 
X(t,e) = {p{t,e),q{t,e)) where q(t,e) crosses exactly m times through the 
horizontal line q = in the interval [0,nT/2]. 
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Proof. Let X(t; ^, 0) = ^, 0), q(t, ^, 0)) a solution (49) in the autonomous 
case (e = 0) satisfying the initial condition 



p(o,e,o) = o, g(o,e,o) = e, 

and the boundary condition 

p{nT/2,^,0) = 0. 



(52) 



(53) 



In Appendix 2 we prove that the family of periodic solutions X{t] ^, 0) of 



(49) for < ^ < /i/cr has an increasing period function T(^) such that 
limT(^) = Tip and lim T(^) = oo, 

with Tip = 271 /u. Moreover, in the autonomous case we have the following 
symmetry 

{t,P,q) {-t,p, -g), 

therefore we can assume that ^ > 0. Given n G N, X(t; ^, 0) is a nT-periodic 
solution of (49) if and only if there is an integer m > 1 such that 

mT{S,) = nT. 

Since 



(54) 



miT{i)=Ti 



ip 



2tt 



CO 



nT 2tt unT 

we have — > — therefore m < . 

m UJ 27r 



Let 



/i/a > ^1 > 6 > ■ ■ • > > 0, 



be the solutions of (54) with m = 1,2, 
dary condition (53) Then 



,'dn- Since we consider the boun- 



/Co = {eGM:p(nT/2,e,0) = 0} 



Now we compute the index for G /Co- Note the following 

T(o ^ T(6: 



If ^ < a then 



< 



: nT/2, 
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. Ife>ei then ^ > ^ = nT/2. 

Since T{C,) is increasing, we have for values ^ close to .^i that 

p(nT/2;e,0) >0 if ^<a, 
p(nr/2;e,0) <0 if ^ > ^i- 

From here, the Brouwer index ind(p(nT/2, -, 0), .^i) = —1 (see for instance 
|34j for definition and basic properties). From the previous calculus we can 
conclude that in general 

ind((p(nT/2,-,0),e^) = (-l)™. (55) 

By symmetry the indices of —^m are ind(p(nT/2, -, 0), ^m) = (— 1)™"^^ We 
also calculate the index at = 0. We do this by linearization, i.e., 

ind((p(nT/2, -, 0), 0) = sign(|^(nT/2; 0, 0)) . 

To this end we consider the linearized problem of at (0, ^) and we observe 
dX 

that -g^(t,^,0) is solution of the initial value problem 



X = AX, X(0) = (^^ 



therefore 



ind((p(nT/2,-,0),0) =sign(^sin(^)) = (-1)* 

In conclusion for e = there exists -^n nontrivial riT-periodic solutions of 
(g with 

p(0;e,0) = 0, g(0;e,0)=e>0, 

such solutions can be labeled according to the number of times m that the 
function g(t;^,0) passes through the horizontal line g = in [0,nT/2] with 
m = 1, . . . and initial conditions 

/i/a > gi(0) = ^1 > ■ ■ • > g^JO) = > 0. 
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From (55 ) and the Implicit Function theorem, for each n > n^: with n^: 



the is a function E : [0, 
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e -E'(e) such that X{t, e) is a nT- 



periodic solution of (49) that satisfy the initial condition 

X(0,e) = 





m 

and the boundary condition 

p{nT/2,E{e),e) = 0, 
for all e G [0, e„_m[- This completes the proof. 



□ 



6 Conclusions and outlook 

In this work we have shown the appearance of chaotic and periodic behavior 
in chemical systems as a direct consequence of internal fluctuations. We have 
concentrated on the simple reaction A < — > 2A, which at the mean-field (fluc- 
tuations free) level simply shows logistic growth. While this theory correctly 
predicts the short time behavior, long times are dominated by fluctuations. 
In this case we have seen that fluctuations may sustain metastable states and 
periodic orbits. For this same reaction, if we allow the reaction rates to vary 
periodically in time we find the presence of chaotic orbits sustained by chem- 
ical noise, while the mean-field theory only reflects periodic orbits. We have 
also been able to rigorously prove the existence of even and periodic solutions 
of the nonautonomous system as a result of the global continuation of even 
and periodic solutions of the autonomous system. It is important to remark 
here that the deterministic trajectories studied here will not be observed 
in experiments, but their noisy counterparts. That is, small Gaussian dis- 
tributed fluctuations about the deterministic trajectories studied herein will 
be present at any time. We also note that the appearance of periodicity and 
chaos in chemical reactions due to intrinsic fluctuations has already been in- 
vestigated |35l |36l |37] . Nevertheless our approach is fundamentally different 
as these previous studies focused on systems close to a bifurcation threshold. 
So the periodic and chaotic behaviors were already present in the mean-field 
deterministic dynamics, and the internal noise anticipated the threshold. In 
our case, the mean-field dynamics were of too low dimensionality for showing 
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periodic and chaotic orbits respectively. In this sense, the oscillations and 
chaos were purely sustained by chemical fluctuations. 




Figure 6: Optimal paths to extinction in Lotka-Volterra dynamics (57) for 
/i = (7 = A = 1. The initial numbers of predators (represented by the red 
upper line) and preys (represented by the green lower line) are respectively 
NAit = -1) ^ 0.39 and Nsit = -1) ^ 0.15. The extinction takes place at 
t = and for Pa = Pb = 0, Qa = —I and qb = 1. The trajectories lie in if = 1 
manifold. 



We have also studied the optimal paths to extinction in a plankton pop- 
ulation dynamics model. Although extinction phenomena have been consid- 
ered numerous times within this framework [131 El EO], the approach was 
based on the Hamiltonian dynamics on the stationary H = manifold. The 
Hamiltonian dynamical system in our case was however degenerated and all 
the extinction trajectories fell on the H manifolds. We leave for fu- 
ture work the extension of our results to multispecies reactions, which are 
characterized by large deviation Hamiltonian dynamical systems of higher 
dimension. Notably, there is an important problem of a two species reaction 
which is very much related to the plankton extinction described herein. It 
is the appearance of extinction events in the predator-prey Lotka-Volterra 
dynamics. This dynamics is formalized by means of the reaction set 

B^2B, A + B^2A, (56) 

where A is the predator and B is the prey. This system has been stud- 
ied by means of a diffusion approximation and solving the corresponding 



26 



Fokker- Planck equation [3H]- A different alternative is using a large devia- 
tions approach, which yields the Hamiltonian 



^ = KPa - l)ga - O-pbiPb - l)gb + Apa(j9fe - Pa)qaqb- (57) 

The numbers of predator and prey are given respectively by NA{t) = Pa{t)qa{t) 
and NB(t) = Pbit)qb(t). The optimal paths to extinction appear again on 
H ^ manifolds. We have numerically computed two of them in Fig. |6| A 
large deviation drives the system to extinction which takes place when pa 
and pb become zero. One can see in this figure that the decreasing number 
of preys pulls the predators to extinction, as it is reasonable to expect. Al- 
ternatively, a large deviation driving the predators to extinction will lead to 
an unbounded growth of the preys. We expect that a systematic study of 



the extinction trajectories of the Lotka-Volterra Hamiltonian (57) will shed a 
valuable light on the probabilistic structure of extinction events in predator- 
prey dynamics. 
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Appendix 1 

The purpose of this appendix is to show more exphcitly the calculations to 



obtain the expression of the operator T in the formula (51). To this end, 
consider the system 

X = AX + B{t,X) (58) 

with 

The fundamental matrix of the associated autonomous system X = AX is 
the exponential matrix e*"^ given by 

cos ut sm ut 

= ( 4a; 

— sin ut cos ut 
a 



Then by the method of variations of constants the general solution of (58) is 
given by 

X{t) = e'^C+ [ e^'-'^^B{s,X)ds, (59) 
Jo 

where C a constant vector. Since we are looking for T-periodic solutions 
imposing the boundary conditions X(0) = X(T) we get 

C = e^^C+ / e'^'-'^^B{s,X)ds 
Jo 

this implies 

{h-e^^)C= r e^''-'^''B{s,X)ds 
Jo 



C=ih- e 
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T 

TA 



1 f e(^-^)^5(s,X)rfs, 
Jo 



replacing this in the formula (59) we obtain 
On the other hand 



Jo Jo 



a 



TA\ 



J ={12- e 



since det J = 2(1 — coswT) then 



1 — cos uT — sm uT 

sm uT 1 — cos uT 

a 



1 



1 

4w cos(wT/2) 



\ a sin(wT/2) 
by direct calculation it is found that 



cr cos(a;T/2)\ 

I 



Auj sin(u;T/2) 



Therefore 



X{t) = J-\^^ r e^''-'^^B{s,X)ds+ t e^'-'^^B{s,X)ds 
Jo Jo 

X{t) = [ J-^e^^ + h] f e(*-^)^5(s, X) ds + J-^e^^ / e(*-'^)^5(s, X) ds 

Jo Jt 

The matrix J satisfies J^^ = J^^e^^ + l2- In consequence 

X{t) = J-' f e^'-'^^B{s,X)ds + J-\^^ r e^'--'^^B{s,X)ds 
Jo Jt 



We define the Green's matrix 



G{t,s) -- 
Thus we can write 



j-l^TA^it-s)A. 0<t<S<T 



X{t) = [ G{t,s)B{s,X)ds 
Jo 
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Let Gi — J ^e^^ Explicitly this matrix is given by 



cosLc{l-s) siii(^:/72)-c()s(^:/72) sin ■j(f-s) 

2sin(a;T/2) J 



CT rsina,'(/ — s) siu(^'772)+cxjs(^':/72) cos ,s) 
"iJ ^ 2sin(^:/72j 

COS^'l^/ — S'j siii(^''i72)— cus(i^''i72) hU\U)(t — s) 
2 sin(£jr/2) 



Prom the trigonometric identities follows easily 



/ sin^'(772 - (f - s)) 

sin(a;T/2) 
4a;cosa;(T/2 - {t - s)) 

\a sin(a;T/2) 



0" cos 



sc.(7y2-(f-,s))\ 



4a; sin(cjT/2) 
sinw(T/2 - {t- s)) 



sin(a;T/2) 



Let G2 — J ^e^^e^* Explicitly this matrix is given by 

1 a cos(tjT/2) 



4c.2sin(cur/2) n cosw(t-s) -^sina;(t-s) 
2 ^ ^ 1^ sinc(;(t — s) cosa;(i — s) 



^2 — 1 cos(ajT/2) 1 | I 4<^ 

(T 2sin(a;T/2) 

By direct calculation 

- cosa)(t— s) sin(tjT/2)— cos(t<jT/2) sinu)(t— s) cr rcos(aiT/2) cosu)(t— s)— sina;(t— s) sin(u;r/2) 

/-( _ I 2sin(wT/2) ~4u;l. 2sin(tjT/2) 

^-^2 — I r cos{ojT/2) cosuj{t-s)-smuj(t-s) sin(u;T/2) ] - cosa;(t-s) sin(t^T/2)-cos(a;T/2) smoj{t-s) 



2 sin(a;r/2) 



2 sin(a;r/2) 



again, using basic trigonometric identities follows that 

/ sina;((s-t) -T/2) cosa;((s - t) - T/2) \ 



J-l^TA^{t-s)A 



sm{ujT/2) 4uj sin(wT/2) 

4a;cosa;((s-t) - T/2) sina;((s-t) - T/2) 



\ a sm{ujT/2) 



sm{ujT/2) 



I 



In orther to present a upper bound for the constant P{uj,a) we give the 
calculations for 



max / \Gii(t,s)\ds 



To this end, observe that 



^0 



\Gi^i{t,s)\ds 



1 



^0 



sina;(T/2 - {t - s)) 



sm{uT/2) 



ds- 



smuj{T/2- (s-t)) 



s\n{uT/2) 



ds 
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On the other hand 

1 /•* sinw(T/2 - (t-s)) 







sin(u;T/2) 



ds 



Juj(T/2-t) 



smu 



sm{uT/2) 



and also 

sina;(r/2- (s - t)) 



sin(wT/2) 







sinw 


(is = 




sin(wr/2) 


2u 


luiT/2 



du. 



It follows that 



a;r/2 



sm-u 



(is 



,0 ^ Jo |sin(u;T/2)| 

The same analysis applies to other components of the matrix G and we obtain 

r-T 







|G'i,2(t, s)\ds 



\G2,i{t, s)\ds 



|G'2,2(i, s)|cis 







sin(a;r/2)| 



4 I cos u I 

(J Jo I sin(a;T/2)| 

]^ /•ajT/2 



smw 



ds. 



Since 



|sin(a;T/2)| 

B = B(a,a)= max / \G(t,s)\ds= max ( max / \Gi j(t, s)\ds] . 
follows that 

Appendix 2 

In this appendix we are going to study the period function of periodic so- 
lutions around the center {1/2, fi/a) for the autonomous case. Given the 



T 



—— max > , , 

sina;r/2| \ 2 a 8u 



1 2 a 



clear symmetry of the vector field of our Hamiltonian system (37) around 
the center we consider the next change of variables 



(60) 
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We obtain in the new coordinates (p, q) the quadratic Hamiltonian function 



cr 



(t2 



(61) 



^{Piq) has the following symmetries 
The corresponding dynamical system is 



Si 
S2 
S3 



cr 



-q + crgp 



^ - 2 

— p — apq 

a 



(62) 



The invariant lines are now 

{p=-l/2}, {p = l/2}, {q 



-/i/cr} and {? = yu/cr} 



They define a new quadrangular area A which center is our equilibrium point 
(0,0). See Fig.|7l 




Figure 7: Quadrangular region A. 



Consider a periodic solution X{t) = {p{t),q{t)) of (62) inside A with 
initial conditions 

P(0) = C, g(0) = 0, (63) 
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this implies that ( satisfies —l/2<(< 1/2. Since "H is a first integral of 



(62) we take the level set 



and consider the equation 



which can be written as 



'Hm,q{t)) = h 



1, 



By the symmetry S2 we can take < < < 1/4 and get from direct 
calculus 



ant) 



4p2(t) _ 1 ' 



From (62) it follows that 



P 



and we finally have 



p\t) = \^^\p\t)-e){4p\t)-l). 

Note that the right hand side in the last equation is positive. Let T{() be 
the period of the solution X(t) of (62) satisfying (63). We have 



T{C)/4 



dt 



no 



U p = (v, then 



T(C)/4 



dp 



dt. 



T(C) = -K{20, 



dp. 



35 



where K{x) 



1 



V(l -xV) 



dv is the complete elhptic integral of 



the first kind. On the other hand the linearized problem of (62) about the 
equilibrium solution (0, 0) is given by 



a 



a 



1/2 



yi 



and its general solution is 

^, = -C " sin(^t)+C,cos(|t 



7/2 = Ci COS (J^t^ + sin (^f^t 



Air 



This shows that the period of the linear problem is Tip 

From the above discussions we are now in position to prove the next 
statement. 



Proposition 1. Let X{t;() = {p{t;C),q{t]()) be a non trivial periodic so- 
lution of (62) that satisfies the initial conditions (63). IfT{Q is the period 
function of X{t; () this function satisfies the following properties 

(a) ^hm T(C)=T,,. 

(b) lim T(C) = oo. 
dT 

(c) — > for allO<C <k- 
dC 2 



Proof. The proof of (a) and (6) follows from the well know properties of the 
complete elliptic integral of the first kind function. Moreover 

K(2n = / dv = / =de. 
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To prove (c) consider the real function 

mc) = 



1 - 4^2 gin2 Q 

with / : [0,7r/2] x [0,l/2[^ M. This function satisfies 

• /(■ , C) is a Riemann integrable function in [0, vr/2], for all G [0, l/2[, 

• f{6, ■ ) is a different iable function in [0, l/2[, for all 6 G [0, vr/2], 

df 

furthermore — is continuous in [0, 7r/2] x [0, l/2[. The previous observations 

over the function / are the hypothesis of the classic version of the rule of 
derivation under the integral, then 

^ ' > 0, 



dC fx Jo (1 -4C2sin'^)3/2 
and this proves (c). □ 



Observation. By symmetry of the system (62) , the above proposition is 
also true if we consider non trivial solutions Xit^rj) = {p{t;r)),q{t;rj)) that 
satisfies the initial conditions 

p(0) = 0, g(0) = T], 

with < ?7 < /i/cr. 
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